home *** CD-ROM | disk | FTP | other *** search
/ MacWorld 1997 September / Macworld (1997-09).dmg / Serious Software / Cherwell Scientific Demos / pro Fit / pro Fit 5.0 demo (fpu).sea / pro Fit 5.0 demo (fpu) / Functions & Programs / •Gadgets / Projection 3D < prev    next >
Text File  |  1996-06-01  |  4KB  |  117 lines
  1. {
  2.   this program takes a three dimensional data set, defined by its x-, y-, and z- coordinates, and projects
  3.   its points onto the plane perpendicular to a vector e  in such a way that the position vector of 
  4.   the first point (i.e. the one in the first non empty row of the x-y-z columns) that is read in is 
  5.   made vertical. 
  6.   The two dimensional picture obtained by the projection is then rotated by an angle phi.
  7.   
  8.   The resulting data set is stored in two data columns.
  9.   
  10.   To make things faster, the program assumes that if the xCol contains a valid
  11.   number, also yCol and zCol do.
  12.   
  13.   To use this program, choose "Add To Menu" from the Misc menu (or click the Add button).
  14.   Then choose the program's name from the Misc menu.
  15.  
  16. }
  17.  
  18. program Projection3D;
  19.  
  20. var e1,e2,e3,phi,xCol,yCol,zCol:extended;
  21.             a1,a2,a3,b1,b2,b3,c1,c2,c3:extended;
  22.             xx,yy,xxC,yyC,L,phi0, CosPhi, SinPhi:extended;
  23.             i:integer;
  24.  
  25. procedure initialize;
  26. begin
  27.   e1:=1;e2:=1;e3:=3;phi:=0;
  28.   xCol:=1;yCol:=2;zCol:=3;xxC:=4;yyC:=5;
  29. end;
  30.  
  31. function sp(x1,x2,x3,y1,y2,y3:extended):extended;    {scalar product of two vectors}
  32. begin
  33. sp:=x1*y1+x2*y2+x3*y3
  34. end;
  35.  
  36. begin
  37. {reset phi. See below}
  38. phi:=(phi-phi0)*180/pi;
  39. {first, ask for input and output columns}
  40. SetBoxTitle('Input and output columns');
  41. for i:=1 to NrCols do 
  42. begin
  43.     if ColEmpty(i) then leave
  44. end;
  45. xxC:=i;
  46. for i:=xxC+1 to NrCols do 
  47. begin
  48.     if ColEmpty(i) then leave
  49. end;
  50. yyC:=i;
  51. Input('$Cx Column',xCol,'$Cy Colum',yCol,'$Cz Colum',zCol,'$Coutput x Col',xxC,'$Coutput y-col',yyC);
  52.  
  53. {then, ask the viewing direction (e) and the rotation angle of the final 2-D output}
  54. SetBoxTitle('projection direction and rotation angle');
  55. Input('e1=',e1,'e2=',e2,'e3=',e3,'rotation angle',phi);
  56.  
  57. {normalise the vector e}
  58.         L:=sqrt(sqr(e1)+sqr(e2)+sqr(e3));
  59.         e1:=e1/L;
  60.         e2:=e2/L;
  61.         e3:=e3/L;
  62.  
  63. {find a unit vector a perpendicular to e}
  64. if e2=0 then begin a1:=0;a2:=1;a3:=0 end else
  65.     begin
  66.             a1:=1;a2:=-e1/e2;a3:=0;
  67.             a1:=a1/sqrt(1+sqr(a2));
  68.             a2:=a2/sqrt(1+sqr(a2));
  69. end;
  70.  
  71. {calculate the unit vector perpendicular to both e and a: cross product!}
  72. b1:=e2*a3-e3*a2;
  73. b2:=e3*a1-e1*a3;
  74. b3:=e1*a2-e2*a1;
  75.  
  76. i:=0;
  77. repeat i:=i+1 until dataok(i,xCol); {find the first data row}
  78.  
  79. {find the length L of the projection of the position vector on e}
  80.     L:=sp(data[i,xCol],data[i,yCol],data[i,zCol],e1,e2,e3);
  81.     
  82. {find c, the component of the position vector perpendicular to e}
  83.     c1:=data[i,xCol]-L*e1;c2:=data[i,yCol]-L*e2;c3:=data[i,zCol]-L*e3;
  84.     
  85. {find the lenghts xx and yy of the components of c along the unit vectors a and b}
  86.     xx:=sp(c1, c2,c3,a1,a2,a3);
  87.     yy:=sp(c1, c2,c3,b1,b2,b3);
  88.     
  89. {the first position vector read in is made vertical by default, do not do it if the projection of the position vector gives zero}
  90.     if (xx=0) and (yy=0) then phi0:=0 else phi0:=pi/2-arctan(yy/xx);
  91.     {rotating th output by phi0 produces verticality of the first position vector which is read in}
  92.     {add the user defined angle to this angle}
  93.     phi:=pi*phi/180+phi0;
  94.     
  95.     {calculate the elements of the rotation matrix}
  96.     CosPhi:=cos(phi);SinPhi:=Sin(phi);
  97.     
  98. {store the rotated output in the two output columns}
  99.     data[i,xxC]:=xx*CosPhi-yy*SinPhi;
  100.     data[i,yyC]:=xx*SinPhi+yy*CosPhi;
  101.     
  102.     i:=i+1;{go to the next row }
  103.  
  104. while (i<=NrRows)  do    {calculate the projection for each data raw}
  105. begin
  106. if  dataok(i,xCol) then
  107. begin
  108.     L:=sp(data[i,xCol],data[i,yCol],data[i,zCol],e1,e2,e3);
  109.     c1:=data[i,xCol]-L*e1;c2:=data[i,yCol]-L*e2;c3:=data[i,zCol]-L*e3;
  110.     xx:=sp(c1, c2,c3,a1,a2,a3);
  111.     yy:=sp(c1, c2,c3,b1,b2,b3);
  112.     data[i,xxC]:=xx*CosPhi-yy*SinPhi;
  113.     data[i,yyC]:=xx*SinPhi+yy*CosPhi;
  114. end;
  115.     i:=i+1;
  116. end;
  117. end;